ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.
Sujet2-Etape1
1 Etape 1 : Première manipulation du RPG
1.1 Je choisis un point sur la carte : https://www.google.fr/maps
=> clic droit sur la carte + clic gauche sur les coordonnées (copiées dans le presse-papier)
J’indique un rayon en mètres (< 15000) pour sélectionner les parcelles autour du point
Code
# coord_gmaps<- rstudioapi::showPrompt(title = "Collez les coordonnées", message = "coordonnées Gmaps", default = "")
coord_gmaps<-"46.46946179131805, -1.4932577154775046"
lat<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,1])
lon<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,2])
# rayon<-as.numeric(rstudioapi::showPrompt(title = "Rayon", message = "Rayon (en m)", default = ""))
rayon<-5000
# création d'une table sf «point» avec les coordonnées saisies
# transformation des coordonnées en syst de proj 2154 (Lambert II - Français)
point<-data.frame(lon,lat,rayon) %>%
st_as_sf(coords = c("lon","lat"),crs = "EPSG:4326") %>%
mutate(coord_pt_gps=st_as_text(geometry)) %>%
st_transform("EPSG:2154") %>%
st_sf() %>% clean_names() %>%
rename(geom=geometry)
# st_crs(point)
# st_geometry_type(point)
# plot(point)La table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL/PostGIS.
1.2 Je me connecte au serveur PostgreSQL avec le mot de passe disponible sur Onyxia/Mes services/PostgreSQL/Read me
Code
# le mot de passe est stocké dans un secret Vault
# postgresql_password <- rstudioapi::askForPassword(prompt = "Entrez le password PostgreSQL")
postgresql_password <- "1tfawt3nj7fgzo3w7cma"
# Connection à PostgreSQL
cnx <- dbConnect(PostgreSQL(),
user = "projet-funathon",
password = postgresql_password,
host = "postgresql-438832",
dbname = "defaultdb",
port = 5432,
options="-c search_path=rpg,public") # specify what schema to connect to1.3 je lance une requête SQL pour sélectionner les parcelles situées autour de mon point dans le rayon choisi.
Code
# suppression de la table «point» si elle existe
dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.point CASCADE;")<PostgreSQLResult>
Code
# écriture de la table point dans la base PostGis
write_sf(point, cnx, append = T)
# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.point ADD CONSTRAINT point_pkey PRIMARY KEY(coord_pt_gps);")<PostgreSQLResult>
Code
# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")<PostgreSQLResult>
Code
# dbExecute(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")
# 3 - Exécution de la requête de découpage du RPG autour du point sur PostGis -------
dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.parc_prox CASCADE;")<PostgreSQLResult>
Code
dbSendQuery(cnx,"CREATE TABLE rpg.parc_prox AS
SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.*
FROM rpg.point p, rpg.parcelles r
WHERE ST_DWithin(p.geom,r.geom,p.rayon);")<PostgreSQLResult>
Code
# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.parc_prox ADD CONSTRAINT parc_prox_pk PRIMARY KEY(id_parcel);")<PostgreSQLResult>
Code
# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX parc_prox_geom_idx ON rpg.parc_prox USING gist(geom);")<PostgreSQLResult>
Code
# 4 - Téléchargement des parcelles proches sous R-------------------------------
# parc_prox<-read_sf(cnx,parc_prox)
parc_prox<-st_read(cnx, query="select * from rpg.parc_prox;")1.4 J’affiche les parcelles autour de mon point avec une carte interactive leaflet
Code
#plot(st_geometry(parc_prox))
#plot(st_geometry(point),add = T,col = "red")
# Transformation de la projection car leaflet ne connait que le WGS 84
carte_parc_prox<- parc_prox %>% st_transform(4326)
# Marqueur du point
pt_mark<- point %>% st_transform(4326)
# ajout du libellé des cultures
carte_parc_prox_lib<-carte_parc_prox %>%
left_join(lib_cult %>% select(-code_groupe_culture),by=c("code_cultu"="code_culture"))
#création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la carte.
labels<-sprintf("<strong>id_parcel : </strong>%s<br/>
<strong>Groupe culture : </strong>%s<br/>
<strong>Culture : </strong>%s<br/>
<strong>Surface (ha) : </strong>%s<br/>
<strong>Département : </strong>%s<br/>
<strong>Commune : </strong>%s<br/>",
parc_prox$id_parcel,carte_parc_prox_lib$libelle_groupe_culture,
carte_parc_prox_lib$libelle_culture,parc_prox$surf_parc,
parc_prox$insee_dep,parc_prox$nom_com) %>%
lapply(htmltools::HTML)
# labels
# création d'une palette de couleurs associée au groupe de culture
factpal <- colorFactor("Paired", parc_prox$code_group)
carte_parc_prox_html <- leaflet(carte_parc_prox_lib) %>%
# addProviderTiles("Esri.WorldImagery") %>%
# addTiles() %>%
addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%
addPolygons( #fillColor="white",
fillColor=~factpal(code_group),
weight = 2,
opacity = 1,
color = "#ffd966",
dashArray = "3",
fillOpacity = 0.5,
highlight = highlightOptions(
weight = 5,
color = "#A40000",
dashArray = "",
fillOpacity = 0.0,
bringToFront = TRUE),
label = labels,
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto",
encoding="UTF-8")) %>%
addMarkers(data=pt_mark,~lon,~lat, popup = ~coord_pt_gps, label = ~coord_pt_gps)Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Paired is 12
Returning the palette you asked for with that many colors
Code
# leaflet(pt_mark) %>% addTiles() %>% addMarkers(~lon,~lat)
carte_parc_prox_htmlCode
# saveWidget(widget = carte_parc_prox_html, file = "carte_parc_prox_html.html")1.5 Quelle est la structure des parcelles agricoles autour de mon point ?
Code
t1 <- parc_prox %>% st_drop_geometry() %>% count(code_group) %>%
add_tally(n) %>%
mutate(n_pct=round(100*n/nn,1)) %>%
select(-nn) %>% rename(n_parcelles=n) %>%
# adorn_totals() %>%
cbind(
# comptage des surfaces
parc_prox %>% st_drop_geometry() %>% count(code_group,wt=surf_parc) %>%
add_tally(n) %>%
mutate(surf_pct=round(100*n/nn,1)) %>%
select(-nn) %>%
rename(surf_parc_ha=n) %>% select(surf_parc_ha, surf_pct) # %>%
# adorn_totals()
) %>% left_join(lib_group_cult,by=c("code_group"="code_groupe_culture")) %>%
select(code_group,libelle_groupe_culture,everything()) %>%
# arrange(as.numeric(code_group)) %>%
arrange(desc(surf_parc_ha)) %>%
adorn_totals() %>%
mutate(taille_moy_parc=round(surf_parc_ha/n_parcelles,1))Storing counts in `nn`, as `n` already present in input
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.
Code
t1 %>%
setNames(c("code","groupe de cultures","nombre de parcelles","(%)","surface (ha)","surface (%)","taille moyenne (ha)")) %>%
kable(
format="html",
caption="<span style='font-size:medium'>Groupes de cultures <strong>locales</strong> par surfaces décroissantes</span>",
format.args = list(decimal.mark = ",", big.mark = " "),
booktabs = TRUE) %>%
kable_styling(font_size = 15) %>%
gsub("font-size: initial !important;",
"font-size: 20pt !important;",.)%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
row_spec(nrow(t1), bold = T, color = "white", background = "grey")| code | groupe de cultures | nombre de parcelles | (%) | surface (ha) | surface (%) | taille moyenne (ha) |
|---|---|---|---|---|---|---|
| 2 | Maïs grain et ensilage | 306 | 18,8 | 1 669,06 | 30,4 | 5,5 |
| 18 | Prairies permanentes | 402 | 24,7 | 1 168,78 | 21,3 | 2,9 |
| 1 | Blé tendre | 178 | 10,9 | 737,43 | 13,4 | 4,1 |
| 19 | Prairies temporaires | 175 | 10,7 | 501,54 | 9,1 | 2,9 |
| 4 | Autres céréales | 88 | 5,4 | 428,99 | 7,8 | 4,9 |
| 6 | Tournesol | 76 | 4,7 | 319,70 | 5,8 | 4,2 |
| 16 | Fourrage | 39 | 2,4 | 141,23 | 2,6 | 3,6 |
| 5 | Colza | 31 | 1,9 | 133,49 | 2,4 | 4,3 |
| 7 | Autres oléagineux | 13 | 0,8 | 74,82 | 1,4 | 5,8 |
| 11 | Gel (surfaces gelées sans production) | 86 | 5,3 | 73,82 | 1,3 | 0,9 |
| 25 | Légumes ou fleurs | 22 | 1,4 | 70,15 | 1,3 | 3,2 |
| 28 | Divers | 180 | 11,1 | 65,18 | 1,2 | 0,4 |
| 3 | Orge | 20 | 1,2 | 52,31 | 1,0 | 2,6 |
| 15 | Légumineuses à grains | 9 | 0,6 | 30,65 | 0,6 | 3,4 |
| 8 | Protéagineux | 3 | 0,2 | 19,63 | 0,4 | 6,5 |
| Total | - | 1 628 | 100,1 | 5 486,78 | 100,0 | 3,4 |
Code
# rm(t1)1.6 Comparaison avec la répartition des cultures au niveau départemental et national
Code
# jointure spatiale avec la couche département pour récupérer le département où tombe le point
# ouvrir la couche des communes (à convertir en Lambert II epsg 2154)
dep <- s3read_using(
FUN = sf::read_sf,
layer = "departement",
object = "2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",
bucket = "projet-funathon",
opts = list("region" = "")) %>%
st_transform(2154)
# jointure
df<-point %>% st_join(dep) %>% st_drop_geometry() %>% select(insee_dep)
dep_pt<-df[1,1]
rm(df)
# sinon sélection du département regroupant la plus grande surface agricole
# cas où le cercle recouvre plusieurs départements
# df<-parc_prox %>% st_drop_geometry() %>%
# count(insee_dep,wt=surf_parc) %>%
# arrange(desc(n))
# dep_pt<-df[1,1]
# rm(df)
# calcul des % surfaces autour du point
stat_pt <- parc_prox %>% st_drop_geometry() %>%
count(code_group,wt=surf_parc) %>% add_tally(n) %>%
mutate(pct_surf_local=round(100*n/nn,1)) %>%
select(code_group, pct_surf_local) Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.
Code
# récup des % surfaces départementales
stat_dep_pt<-s3read_using(
FUN=readr::read_rds,
object = "2023/sujet2/diffusion/resultats/stat_group_cult_by_dep.rds",
bucket = "projet-funathon",
opts = list("region" = "")) %>%
filter(insee_dep %in% dep_pt) %>%
select(insee_dep,code_group,libelle_groupe_culture,pct_surf) %>%
rename(pct_surf_dep = pct_surf)
# récup des % surfaces nationales
stat_fm<-s3read_using(
FUN=readr::read_csv,
object = "2023/sujet2/diffusion/resultats/stat_group_cult_fm.csv",
col_types = cols(code_group = col_character()),
bucket = "projet-funathon",
opts = list("region" = "")) %>%
select(code_group,libelle_groupe_culture,pct_surf) %>%
rename(pct_surf_fm = pct_surf)
# appariement des stas locales, départementales, nationales
stat_compar<-stat_fm %>%
left_join(stat_dep_pt %>% select(code_group, pct_surf_dep),by="code_group") %>%
left_join(stat_pt,by="code_group") %>%
select(libelle_groupe_culture,pct_surf_local,pct_surf_dep,pct_surf_fm) %>%
arrange(desc(pct_surf_local)) %>% adorn_totals()
stat_compar %>%
setNames(c("Groupe de cultures","surf. locales (%)","surf. départ. (%)","surf. France m. (%)")) %>% kable(
format="html",
caption="<span style='font-size:medium'>Comparaison des surfaces locales, départemenales et nationales</span>",
format.args = list(decimal.mark = ",", big.mark = " "),
booktabs = TRUE) %>%
kable_styling(font_size = 15) %>%
gsub("font-size: initial !important;",
"font-size: 20pt !important;",.)%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
row_spec(nrow(stat_compar), bold = T, color = "white", background = "grey")| Groupe de cultures | surf. locales (%) | surf. départ. (%) | surf. France m. (%) |
|---|---|---|---|
| Maïs grain et ensilage | 30,4 | 19,4 | 10,0 |
| Prairies permanentes | 21,3 | 29,0 | 27,7 |
| Blé tendre | 13,4 | 16,4 | 17,7 |
| Prairies temporaires | 9,1 | 9,5 | 5,3 |
| Autres céréales | 7,8 | 8,1 | 3,9 |
| Tournesol | 5,8 | 3,7 | 2,5 |
| Fourrage | 2,6 | 4,1 | 3,7 |
| Colza | 2,4 | 2,3 | 3,5 |
| Autres oléagineux | 1,4 | 0,4 | 0,7 |
| Gel (surfaces gelées sans production) | 1,3 | 0,6 | 1,6 |
| Légumes ou fleurs | 1,3 | 1,2 | 1,5 |
| Divers | 1,2 | 1,0 | 1,3 |
| Orge | 1,0 | 2,3 | 6,2 |
| Légumineuses à grains | 0,6 | 0,4 | 0,2 |
| Protéagineux | 0,4 | 0,9 | 1,2 |
| Plantes à fibres | NA | 0,3 | 0,4 |
| Riz | NA | NA | 0,0 |
| Estives et landes | NA | 0,1 | 8,1 |
| Vergers | NA | 0,1 | 0,4 |
| Vignes | NA | 0,2 | 2,2 |
| Fruits à coque | NA | 0,0 | 0,2 |
| Oliviers | NA | NA | 0,0 |
| Autres cultures industrielles | NA | 0,1 | 1,8 |
| Total | 100,0 | 100,1 | 100,1 |
1.7 Graphique de comparaison des cultures au niveau local, départemental et national
Code
# je sélectionne les 10 groupes de cultures les plus répandus au niveau local
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>%
rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)
# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name)
# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)
tab_piv[is.na(tab_piv)] <- 0
# levels(as.factor(tab_piv$secteur))
# tab_piv$secteur<-factor(tab_piv$secteur, levels=c("france","departement","local"))
# tab_piv<-tab_piv %>% arrange(desc(secteur), desc(value))
# ggplot => problème de tri : affichage des prairies avant le maïs ?
p<-ggplot(tab_piv, aes(x = reorder(libelle_groupe_culture,+value),
y = value,
fill=factor(secteur, levels=c("france","departement","local")))) +
geom_col(position = "dodge") +
labs(title="Surfaces comparées des 10 principales cultures locales, en %", x="Culture", y = "%", fill = "Secteur") +
theme_classic()
# je flippe le graphique pour avoir des barres horizontales
p+coord_flip()1.8 Je teste un graphique par secteur avec facet_wrap
Code
# je sélectionne les 10 groupes de cultures les plus répandus au niveau local
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>%
rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)
# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name)
# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)
tab_piv[is.na(tab_piv)] <- 0
# ggplot => problème de tri : affichage des prairies avant le maïs ?
ggplot(tab_piv, aes(x = reorder(libelle_groupe_culture,+value),
y = value)) +
geom_col(fill = "lightblue", colour = "black",position = "dodge") +
labs(title="Surface par culture", x="Culture", y = "%", fill = "Secteur") +
geom_text(aes(label = value), hjust = -0.3, size = 8/.pt, colour = "black") +
theme_classic() + coord_flip() +
facet_wrap(~secteur,nrow=3,scales='free')